# -*- coding: utf-8 -*-
"""
Created on Tue Oct 13 08:48:48 2020

Helper functions for distance, speed calculations on points data

@author: rklotz
"""

import arcpy 
import pandas
import numpy as np

def calcDist(df_points,spatial_reference,spatial_reference_proj):
    fc_points_ds = "in_memory/fc_points_ds"
    #fc_points_ds_merged = "in_memory/fc_points_ds_merged"
    df_points['pnt_ID'] = df_points.reset_index().index
    
    # convert to feature class
    arcpy.da.NumPyArrayToFeatureClass (df_points[['pnt_ID','X','Y','OID'] ].to_records(index=False), 
                                fc_points_ds , ('X','Y'), spatial_reference=spatial_reference)
    
    # now extract to df with correct projection
    df_points_proj = pandas.DataFrame( arcpy.da.FeatureClassToNumPyArray(
                                    fc_points_ds ,field_names=['pnt_ID','SHAPE@X','SHAPE@Y','OID_'], where_clause="" ,
                                     null_value = -9999, explode_to_points= False,
                                    spatial_reference = spatial_reference_proj ) )
    df_points_proj = df_points_proj.rename(columns={'SHAPE@X':'X','SHAPE@Y':'Y'})
    arcpy.Delete_management(fc_points_ds)
    
    # distance calc
    df_points_proj[['Xprev','Yprev']] = df_points_proj.groupby(['OID_'])[['X','Y']].shift(1)
    df_points_proj['d_km'] =  np.sqrt ( (df_points_proj.X - df_points_proj.Xprev)**2 + (df_points_proj.Y - df_points_proj.Yprev)**2 ) / 1000
        
    df_points = df_points.merge(df_points_proj[['pnt_ID','d_km']] , on='pnt_ID',how='left')

    return df_points
    

def checkECA_NA(df_points_ds,eca_na_fc,spatial_reference):
    # spatial join for within ECA
    fc_points_ds = "in_memory/fc_points_ds"
    fc_points_ds_merged = "in_memory/fc_points_ds_merged"
    df_points_ds['pnt_ID'] = df_points_ds.reset_index().index
    arcpy.da.NumPyArrayToFeatureClass (df_points_ds[['pnt_ID','X','Y'] ].to_records(index=False), 
                                fc_points_ds , ('X','Y'), spatial_reference=spatial_reference)
   
    # creating outfield
    outfield = arcpy.Field()
    outfield.name = "ECAinds"
    outfield.type = "String"
    outfield.length = 100
    outfield.scale = 0

    fms = arcpy.FieldMappings() # creating field mappings
    fms.addTable(fc_points_ds)      # adding all fields from lines

    fm1 = arcpy.FieldMap()          # adding field map for getting port IDs
    fm1.addInputField(eca_na_fc,'OID')

    fm1.outputField = outfield
    fm1.mergeRule = "Join"
    fm1.joinDelimiter = ","

    fms.addFieldMap(fm1) # adding additional outfield

    arcpy.SpatialJoin_analysis(fc_points_ds, eca_na_fc, fc_points_ds_merged, join_operation='JOIN_ONE_TO_ONE' 
                           , field_mapping=fms ,  join_type='KEEP_ALL' , match_option='WITHIN' )
   
    df_points_ds_merged = pandas.DataFrame( arcpy.da.FeatureClassToNumPyArray(
                                fc_points_ds_merged ,field_names=["pnt_ID","ECAinds"], where_clause="" ,
                                 null_value = np.nan, explode_to_points= False,
                                spatial_reference = spatial_reference ) )
    arcpy.Delete_management(fc_points_ds)
    arcpy.Delete_management(fc_points_ds_merged)
   
    #df_points_ds_merged[['in_eca09','in_eca11']] = df_points_ds_merged.ECAinds.str.split(',',expand=True)

    df_points_ds_merged['in_eca09'] = df_points_ds_merged.ECAinds.str.contains('1', na=np.nan, regex=True)*1
    df_points_ds_merged['in_eca11'] = df_points_ds_merged.ECAinds.str.contains('2', na=np.nan, regex=True)*1
    df_points_ds_merged['onland'] = df_points_ds_merged.ECAinds.str.contains('3', na=np.nan, regex=True)*1
    df_points_ds_merged['studyarea'] = df_points_ds_merged.ECAinds.str.contains('4', na=np.nan, regex=True)*1
    df_points_ds_merged['in_naeca'] = df_points_ds_merged.ECAinds.str.contains('5', na=np.nan, regex=True)*1

    df_points_ds = df_points_ds.merge(df_points_ds_merged[['pnt_ID','in_eca09','in_eca11','onland','studyarea','in_naeca']] , on='pnt_ID',how='left')

    return df_points_ds


#def checkECA(X, Y, CAeca_geom, spatial_reference):
#    # create point geometry
#    pnt_geom = arcpy.PointGeometry(arcpy.Point(X,Y), spatial_reference)
#
#    in_eca = pnt_geom.within(CAeca_geom)*1
#
#    return in_eca

def getDist(Xprev, Yprev, X, Y, spatial_reference):
    # create point geometry
    pnt_prev_geom = arcpy.PointGeometry(arcpy.Point(Xprev,Yprev), spatial_reference)
    pnt_geom = arcpy.PointGeometry(arcpy.Point(X,Y), spatial_reference)

    dist_km = pnt_prev_geom.angleAndDistanceTo(pnt_geom,'GEODESIC')[1] / 1000

    return dist_km


# getting bearing between two points
def bearing(lat1,lon1,lat2,lon2) :
    lat1r,lon1r,lat2r,lon2r = map(np.deg2rad, [lat1,lon1,lat2,lon2])
    
    #const y = Math.sin(λ2-λ1) * Math.cos(φ2);
    #const x = Math.cos(φ1)*Math.sin(φ2) -  Math.sin(φ1)*Math.cos(φ2)*Math.cos(λ2-λ1);
    
    y = np.sin( lon2r - lon1r  ) * np.cos( lat2r );
    x = np.cos( lat1r )*np.sin( lat2r ) - np.sin( lat1r )*np.cos( lat2r ) * np.cos( lon2r - lon1r );
    theta = np.arctan2(y, x)
    brng = np.mod( theta*180/np.pi + 360 , 360 )
    #brng = theta*180/np.pi + 360 % 360 

    return brng


# Define a basic Haversine distance formula
#def haversine(lat1, lon1, lat2, lon2):
#    KM = 6371 # radius of earth at equator
#    lat1, lon1, lat2, lon2 = map(np.deg2rad, [lat1, lon1, lat2, lon2])
#    dlat = lat2 - lat1 
#    dlon = lon2 - lon1 
#    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
#    c = 2 * np.arcsin(np.sqrt(a)) 
#    total_miles = KM * c
#    return total_miles
    
#def iterFun(X,Y,Xprev,Yprev,CAeca_geom09,CAeca_geom11 , spatial_reference) :
#    in_eca09 = checkECA(X, Y, CAeca_geom09, spatial_reference)
#    in_eca11 = checkECA(X, Y, CAeca_geom11, spatial_reference)
#    dist = getDist(Xprev, Yprev, X, Y, spatial_reference)
#    
#    return in_eca09 , in_eca11 , dist 